/* ************************************************************************* **
**    OpenVFIFE - Open System for Vector Form Instrinsic                     **
**                Finite Element Method (VFIFE)                              **
**                GinkGo(Tan Biao)                                           **
**                                                                           **
**                                                                           **
** (C) Copyright 2021, The GinkGo(Tan Biao). All Rights Reserved.            **
**                                                                           **
** Commercial use of this program without express permission of              **
** GinkGo(Tan Biao), is strictly prohibited.  See                            **
** file 'COPYRIGHT'  in main directory for information on usage and          **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.                  **
**                                                                           **
** Developed by:                                                             **
**      Tan Biao (ginkgoltd@outlook.com)                                     **
**                                                                           **
** ************************************************************************* */

// $Date: 2020-05-11 $
// Written: Tan Biao
// Revised:
//
// Purpose: This file contains the class definition for materials, including
// base class BaseMaterial (abstract base class), UniaxialPlastic, NdPlastic,
// and several  derived materials - LinearElastic, NonlinearElastic,
// UniBilinearStatic, UniBilinear, UniPower.

// The interface:
//

#include "../headers/material.h"
using namespace std;

/*****************************************************************************/
/*                                  Helper                                   */
/*****************************************************************************/
inline int sign(int x)
{
    return (x > 0) - (x < 0);
}

inline int sign(double x)
{
    return (x > 0) - (x < 0);
}

/*****************************************************************************/
/*                               BaseMaterial                                */
/*****************************************************************************/
BaseMaterial::BaseMaterial(int id)
{
    id_ = id;
    type_ = "BaseMaterial";
    E_ = 0.0;
    G_ = 0.0;
    nu_ = 0.0;
    density_ = 0.0;
    ce_ = 0.0;
    damp_ratio_ = 0.0;
    damp_coef_ = 0.0;
    mu_ = 0.0;
    limited_tensile_strain_ = 9.99e9;
}

BaseMaterial::~BaseMaterial()
{

}

int BaseMaterial::id() const
{
    return id_;
}

string BaseMaterial::type() const
{
    return type_;
}

double BaseMaterial::E() const
{
    return E_;
}

double BaseMaterial::G() const
{
    return G_;
}

double BaseMaterial::nu() const
{
    return nu_;
}

double BaseMaterial::density() const
{
    return density_;
}

double BaseMaterial::ce() const
{
    return ce_;
}

double BaseMaterial::damp_ratio() const
{
    return damp_ratio_;
}

double BaseMaterial::damp_coef() const
{
    return damp_coef_;
}

double BaseMaterial::mu() const
{
    return mu_;
}

double BaseMaterial::yield_stress() const
{
    return yield_stress_;
}

double BaseMaterial::limited_tensile_strain() const
{
    return limited_tensile_strain_;
}

bool BaseMaterial::isFractured(double strain) const
{
    bool flag = false;
    if (strain > limited_tensile_strain_)
    {
        flag = true;
    }
    return flag;
}

void BaseMaterial::setE(double val)
{
    if (val < 0)
    {
        cerr << "Error: E must be positive!" << endl;
    }
    E_ = val;
}

void BaseMaterial::setG(double val)
{
    if (val < 0)
    {
        cerr << "Error: G must be positive!" << endl;
    }
    else
    {
        G_ = val;
    }
}

void BaseMaterial::setNu(double val)
{
    if (val < 0)
    {
        cerr << "Error: nu must be positive!" << endl;
    }
    nu_ = val;
    G_ = E_ / 2.0 / (1 + nu_);
}

void BaseMaterial::setDensity(double val)
{
    if (val < 0)
    {
        cerr << "Error: density must be positve" << endl;
    }
    density_ = val;
}

double BaseMaterial::calcCe() const
{
    // ce_ = sqrt(E_ / density_);
    // return ce_;
    return sqrt(E_ / density_);
}

void BaseMaterial::setDampRatio(double val)
{
    damp_ratio_ = val;
}

void BaseMaterial::setDampCoef(double val)
{
    damp_coef_ = val;
}

void BaseMaterial::setMu(double val)
{
    mu_ = val;
}

void BaseMaterial::setLimitedTensileStrain(double tensile)
{
    if (tensile > 0.0)
    {
        limited_tensile_strain_ = tensile;
    }
    else
    {
        cout << "Waring: Limited tensile strain must be positive! The input ";
        cout << "tensile is " << tensile << ". And the ";
        cout << "limited_tensile_strain_ will be set as 1e9!" << endl;
    }
}

void BaseMaterial::info() const
{
    cout << "Material ID: " << id_ << endl;
    cout << "Material Type: " << type_ << endl;
    cout << "E: " << E_ << endl;
    cout << "G: " << G_ << endl;
    cout << "density: " << density_ << endl;
    cout << "Ultimate tensile strain: " << limited_tensile_strain_ << endl;
}


/*****************************************************************************/
/*                              LinearElastic                                */
/*****************************************************************************/
LinearElastic::LinearElastic(int id): BaseMaterial(id)
{
    type_ = "LinearElastic";
}

LinearElastic::~LinearElastic()
{

}

double LinearElastic::Eq(double ds, double* stress, double* strain, double* ps,
    double* alpha, double* q) const
{
    (*strain) += ds;
    (*stress) += E_ * ds;
    return E_;
}


/*****************************************************************************/
/*                              NonlinearElastic                             */
/*****************************************************************************/
NonlinearElastic::NonlinearElastic(int id): BaseMaterial(id)
{
    type_ = "NonlinearElastic";
}

NonlinearElastic::~NonlinearElastic()
{

}


/*****************************************************************************/
/*                              UniaxialPlastic                              */
/*****************************************************************************/
UniaxialPlastic::UniaxialPlastic(int id): BaseMaterial(id)
{
    type_ = "UniaxialPlastic";
    m_ = 0.0;
    yield_stress_ = 0.0;

}

UniaxialPlastic::~UniaxialPlastic()
{

}

void UniaxialPlastic::setYieldStress(double val)
{
    if (val <= 0)
    {
        cerr << "Error: yield stress must be positive!" << endl;
    }
    yield_stress_ = val;
}

bool UniaxialPlastic::isYield(double stress, const double* alpha,
    const double* k) const
{
    double f = (stress - (1 - m_) * (*alpha))
               * (stress - (1 - m_) * (*alpha))
               - ((1 - m_) * yield_stress_ + m_ * (*k))
               * ((1 - m_) * yield_stress_ + m_ * (*k));
    // cout << "Yield Function: " << f << endl;
    return f >= -1e-9;
}

void UniaxialPlastic::info() const
{
    BaseMaterial::info();
    cout << "yield stress: " << yield_stress_ << endl;
}


/*****************************************************************************/
/*                                 UniIdeal                                  */
/*****************************************************************************/
UniIdeal::UniIdeal(int id, double E, double yield_stress):UniaxialPlastic(id)
{
    if (E < 0)
    {
        cerr << "Error: E must be positive!" << endl;
    }
    if (yield_stress <= 0)
    {
        cerr << "Error: yield stress must be positive!" << endl;
    }
    type_ = "UniIdeal";
    E_ = E;
    Et_ = 0;
    yield_stress_ = yield_stress;
}

UniIdeal::~UniIdeal()
{

}

bool UniIdeal::isYield(double stress, const double* alpha,
    const double* k) const
{
    double f = stress * stress - yield_stress_ * yield_stress_;
    return f >= -1e-10;
}


double UniIdeal::Eq(double ds, double* stress, double* strain, double* ps,
    double* alpha, double* q) const
{
    /*  Return Mapping Algorithm, << Computational Inelasticity>> Simo */
    // step 1: database at last step, ${{\epsilo}_{n}^{p}, {\alpha}_n, q_n}$
    double Eq = E_;

    // step 2: strain field at current step
    (*strain) += ds;

    // step 3: compute elastic trial stress and test for plastic loading
    double stress_trial = E_ * ((*strain) - (*ps));
    double yield_func_trial = fabs(stress_trial) - yield_stress_;

    // step 4: return-mapping
    if (yield_func_trial < 0)
    {
        // not yield
        (*stress) = stress_trial;
    }
    else
    {
        // yield
        double dgamma = yield_func_trial / E_;
        (*stress) = stress_trial - dgamma * E_ * sign(stress_trial);
        (*ps) += dgamma * sign(stress_trial);
        (*alpha) += 0;
        (*q) += 0;
        Eq = Et_;
    }
    return Eq;
}


/*****************************************************************************/
/*                               UniBilinear                                 */
/*****************************************************************************/
UniBilinear::UniBilinear(int id, double E, double Et, double yield_stress,
    double m):UniaxialPlastic(id)
{
    if (Et < 0 | Et > E)
    {
        cerr << "Error: Et must between 0 and E!" << endl;
    }
    if (yield_stress <= 0)
    {
        cerr << "Error: yield stress must be positive!" << endl;
    }
    if (m < 0 | m > 1)
    {
        cerr << "Error: m must between 0 and 1!" << endl;
    }
    type_ = "UniBilinear";
    E_ = E;
    Et_ = Et;
    yield_stress_ = yield_stress;
    m_ = m;
}

UniBilinear::~UniBilinear()
{

}

double UniBilinear::Et() const
{
    return Et_;
}

void UniBilinear::info() const
{
    UniaxialPlastic::info();
    cout << "Et: " << Et_ << endl;
}

double UniBilinear::Eq(double ds, double* stress, double* strain, double* ps,
    double* alpha, double* q) const
{
    /*  Return Mapping Algorithm, << Computational Inelasticity>> Simo
    */
    // step 1: database at last step, ${{\epsilo}_{n}^{p}, {\alpha}_n, q_n}$
    double Eq = E_;

    // step 2: strain field at current step
    (*strain) += ds;

    // step 3: compute elastic trial stress and test for plastic loading
    double stress_trial = E_ * ((*strain) - (*ps));
    double kesi_trial = stress_trial - (1 - m_) * (*alpha);
    double yield_func_trial = fabs(kesi_trial) - (yield_stress_ + m_ * (*q));

    // step 4: return-mapping
    if (yield_func_trial < 0)
    {
        // not yield
        (*stress) = stress_trial;
    }
    else
    {
        // yield
        double H = E_ * Et_ / (E_ - Et_);
        double dgamma = yield_func_trial / (E_ + m_ * H + (1 - m_) * H);
        (*stress) = stress_trial - dgamma * E_ * sign(kesi_trial);
        (*ps) += dgamma * sign(kesi_trial);
        (*alpha) += dgamma * H * sign(kesi_trial);
        (*q) += dgamma * H;
        Eq = Et_;
    }
    return Eq;
}


/*****************************************************************************/
/*                                 UniPower                                  */
/*****************************************************************************/
UniPower::UniPower(int id, double E, double n, double yield_stress, double m):
    UniaxialPlastic(id)
{
    if (E < 0)
    {
        cerr << "Error: E must be positive!" << endl;
    }
    if (yield_stress <= 0)
    {
        cerr << "Error: yield stress must be positive!" << endl;
    }
    if (m < 0 | m > 1)
    {
        cerr << "Error: m must between 0 and 1!" << endl;
    }
    type_ = "UniPower";
    E_ = E;
    n_ = n;
    k_ = yield_stress / pow(yield_stress / E, n);
    yield_stress_ = yield_stress;
    m_ = m;
}

UniPower::~UniPower()
{

}

void UniPower::info() const
{
    UniaxialPlastic::info();
    cout << "k: " << k_ << ", n: " << n_ << endl;
}

double UniPower::Eq(double ds, double* stress, double* strain, double* ps,
    double* alpha, double* q) const
{
    /*  Return Mapping Algorithm, << Computational Inelasticity>> Simo
    */
    // step 1: database at last step, ${{\epsilo}_{n}^{p}, {\alpha}_n, q_n}$
    double Eq = E_;

    // step 2: strain field at current step
    (*strain) += ds;

    // step 3: compute elastic trial stress and test for plastic loading
    double stress_trial = E_ * ((*strain) - (*ps));
    double kesi_trial = stress_trial - (1 - m_) * (*alpha);
    double yield_func_trial = fabs(kesi_trial) - (yield_stress_ + m_ * (*q));


    // step 4: return-mapping
    if (yield_func_trial < 0)
    {
        // not yield
        (*stress) = stress_trial;
    }
    else
    {
        // yield
        double yield_strain = yield_stress_ / E_;
        // if (fabs(kesi_trial) < 0)
        // {
        //     double sigma = (1 - m_) * (*alpha) - (yield_stress_ + m_ * (*q));
        //     yield_strain = pow(sigma / k_, 1 / n_);
        // }
        // else
        // {
        //     double sigma = (1 - m_) * (*alpha) + (yield_stress_ + m_ * (*q));
        //     yield_strain = pow(sigma / k_, 1 / n_);
        // }
        double Ep = k_ * n_ * pow(fabs(*strain - yield_strain), n_-1);
        double Et = E_ * Ep / (E_ + Ep);
        double H = Et;
        double dgamma = yield_func_trial / (E_ + m_ * Ep + (1 - m_) * Ep);
        // (*stress) = stress_trial - dgamma * E_ * sign(kesi_trial);
        // (*stress) += Ep * dgamma * sign(kesi_trial);
        (*stress) += Et * ds;
        // cout << "Ep * dgamma = " << Ep * dgamma * sign(kesi_trial);
        // cout << ", Et * ds = " << Et * ds << ", E * es = ";
        // cout << E_ * (ds - dgamma * sign(kesi_trial)) << endl;
        (*ps) += dgamma * sign(kesi_trial);
        (*alpha) += dgamma * Ep * sign(kesi_trial);
        (*q) += dgamma * Ep;
        Eq = Et;
        // Eq = H;
    }
    return Eq;
}


/*****************************************************************************/
/*                                 NdPlastic                                 */
/*****************************************************************************/
// NdPlastic::NdPlastic(int id): BaseMaterial(id)
// {
//     type_ = "NdPlastic";
//     Et_.setZero();
//     Eq_.setZero();
//     stress_.setZero();
//     strain_.setZero();
//     delta_strain_.setZero();
// }

// NdPlastic::~NdPlastic()
// {

// }

// void NdPlastic::info() const
// {
//     cout << "Material ID: " << id_ << endl;
//     cout << "Material Type: " << type_ << endl;
//     cout << "E: " << E_ << endl;
//     cout << "G: " << G_ << endl;
//     cout << "Et: " << Et_ <<endl;
//     cout << "density: " << density_ << endl;
// }


// int main()
// {
//     double E = 2.06E11;
//     double Et = E / 100;
//     double yield_stress = 235E6;
//     double m = 0;
//     // UniBilinear mat = UniBilinear(1, E, 2.06E9, yield_stress, m);
//     UniIdeal mat = UniIdeal(1, E, yield_stress);
//     mat.setLimitedTensileStrain(0.003);
//     mat.setNu(0.31);
//     mat.setDensity(7850);
//     // mat.info();

//     double ds = 0.00001;
//     double stress = 0, strain = 0, alpha = 0, k = 0, eps = 0;

//     // alpha = yield_stress / (E * Et / (E - Et));
//     // k = yield_stress;

//     for (int i = 0; i < 1000; i++)
//     {
//         // cout << "\nstep " << i << endl;
//         double E = mat.Eq(ds, &stress, &strain, &eps, &alpha, &k);
//         cout << strain << ", " << stress << endl;

//         // if (mat.isFractured(strain))
//         // {
//         //     cout << "Fractured" << endl;
//         //     break;
//         // }
//     }

//     ds = -ds;
//     for (int i = 0; i < 3000; i++)
//     {
//         // cout << "\nstep " << i << endl;
//         double E = mat.Eq(ds, &stress, &strain, &eps, &alpha, &k);
//         cout << strain << ", " << stress << endl;
//     }

//     ds = -ds;
//     for (int i = 0; i < 2000; i++)
//     {
//         // cout << "\nstep " << i << endl;
//         double E = mat.Eq(ds, &stress, &strain, &eps, &alpha, &k);
//         cout << strain << ", " << stress << endl;
//     }
//     return 0;
// }
